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Low-temperature series are calculated for the free energy, magnetisation, susceptibility 
and field-derivatives of the susceptibility in the Ising model on the quasiperiodic Penrose 
lattice. The series are computed to order 20 and estimates of the critical exponents a, /3 
and 7 are obtained from Pade approximants. 



I. INTRODUCTION 

The problem of the relevance of disorder for phase transitions in lattice models of statistical mechanics 
has attracted attention for many years and the discovery of quasicrystals ffl has served to increase 
interest in the physical properties of disordered systems. A fundamental problem in this field is whether 
quasiperiodic order is strong enough to change the critical behaviour of magnetic phase transitions. To 
investigate this problem we consider in this article a classical Ising model defined on an underlying 
quasiperiodic lattice. 

There have been many works in this field since the late eighties. A heuristic criterion (Harris-Luck 
criterion) has been formulated ^| which relates the critical behaviour to fluctuations of the number of spin 
couplings in a given region. The spatial scaling of fluctuations was described in terms of a "wandering 
exponent" lo which was required to exceed a threshold lo c in order to produce a new universality class. For 
the majority of quasiperiodic structures existing in reality, such as the structural models of quasicrystallinc 
phases discovered so far to can be calculated exactly, due to the self-similarity or inflational symmetry 
of the structure, yielding a value smaller than the threshold and suggesting the irrelevance of disorder. 
However since numerous structures like the rhombic sevenfold or ninefold lattices |3|,^| exist, which are 
deprived of inflational symmetry and are therefore potential candidates for novel critical behaviour, there 
is still a strong motivation for dealing with quasiperiodic Ising models. 

Quasiperiodic Ising models were investigated by Monte-Carlo simulations which at present, seem 
to yield the most precise estimates for the transition temperature and critical exponents. Indeed, in 
Q computations for large periodic approximants (PA) of the Penrose tiling (PT) |9) were carried out 
and obtained values for the correlation length v and the two-spin correlation function r\ exponents with 
two-digit precision (y = 1.02 ± 0.02, 77 — 0.252 ± 0.003) which agreed with the square lattice values 
[y = 1, rj = 0.25). Moreover, the non-universal critical temperature T c has also been determined with an 
impressively small error kT c = 2.398 ± 0.003. 

It is worth mentioning that a novel invaded-cluster algorithm, which modifies the temperature during 
the simulation towards the critical one, as opposed to standard Monte-Carlo algorithms with fixed tem- 
perature, was also applied to quasiperiodic systems ||] to give an improved estimate of T c . The critical 
exponents are not available in this case, however. 

Another approach is an approximate renormalisation group analysis |lO|Jll[| which yields poor results, 
however. For the PT the specific heat exponent equals a — —0.1083 versus a = for the square lattice. 
Quasiperiodic Ising models were also examined by graphical expansion methods |12 L| and by cal- 



culating exact partition functions for PA, obtained from the Kac-Ward determinant |140. In the first 
case estimates for T c and critical exponents have not been considerably improved but this approach 
demonstrated a new feature, a very slow convergence of the partition function (Z) series to its predicted 
asymptotic form. We also investigated the set of zeros of Z in the complex plane (Fisher zeros), which 
turned out to be much more complicated than in the square lattice case. 

The Kac-Ward determinant method appeared to yield highly accurate estimates of the critical temper- 
ature of quasiperiodic Ising models (for example kT c — 2.397820(7) for the PT). Moreover, within this 
framework it was possible to construct a two dimensional Ising model with relevant fluctuations, i.e. for 
which lo > lo c , which shows up another novel feature, namely the divergence of high temperature series 
p6[ . This example is interesting because it points out that in some cases the reliability of methods for 
extracting critical values from analysis of a series expansion, like the Pade- or differential-approximants 
methods, |l5| ] can be questioned. An inspection of the Fisher zeros furnished the explanation, since it 
appeared that the moduli of some complex zeros were smaller than the modulus of the physical singularity 
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(real zero) and thus the complex zeros, rather than the real zero, were limiting the region of convergence 
the series. 

In this paper it is not our purpose to improve on estimating T c or a for quasiperiodic Ising models 
since, due to the slow convergence of series expansions JT4j ] a large, inaccessible number of terms is needed 
to make progress in this field. Instead, we aim at generalising the series expansion approach to the case of 
non-zero field quasiperiodic Ising or Potts models [|l7| and provide alternative estimates of the magnetic 
exponents, /3, and 7. Moreover, this approach allows us to investigate the problem of a disorder-driven 
"softening" of the first-order phase transition in Q-state Potts models for Q > 4 18 1. 



II. THE FINITE LATTICE EXPANSION METHOD FOR ISING MODELS 

The problem consists in calculating the partition function Z(Q) of an Ising model on a lattice Q by 
series expansion. The partition function with field B and coupling constant J is defined in the usual way: 

N , 

Z(G) = ^exp/3{-J^A(a„a fc )-i?^A( ( T„0)} where A( CTl ,a 2 ) = ° (1) 

Wi} U,k) j=i L 

where the sum over spin configurations {dj} = {a±, 02, . . . , &n} consists of N sums each of which runs over 
(Tj = {1,2}. Starting from cluster integral theory (page 42-46 and page 73 in one can formulate a free 
energy (F) expansion in terms of connected graphs for a wide range of models from statistical mechanics. 
In particular, for the non-zero field Ising model or the Q-state Potts models (the generalisation of the 
former one with Q values of spin at each site) the expansion on a lattice Q reads 

\ogZ(G) = ^2(C r ;G)k r {w) where w = tanh/3J (2) 

T" 

where the sum on the right-hand side runs over connected graphs C r from Q. The quantity (C r ; Q) denotes 
the embedding number of C r in Q, counting the number of can be embedded in Q. Finally, the 

weight functions k r {w) depend only on C r not on G. Making use of the independence of weights k r {w) 
from the lattice we can write equation (|^) substituting each connected graph C r for Q, solve the system 
of equations for the weights and plug in the results to equation (|2|). We obtain 

log Z(G) = Y,a r \ogZ(g r ) (3) 

where a r = ^2 p (g P ' 7 G)b rtP and b r ^ p is inverse to the matrix of embedding numbers, i.e. b TiP = (g r ; g p ) ■ 
The sum on the right-hand side in (^) runs over g r from a subset of all connected graphs. It turns out p(| 
that the graph g r can furnish a non- vanishing contribution, i.e. a r ^ 0, if and only if it is an overlap of the 
embeddings of two other graphs having non- vanishing contributions. The construction of graphs therefore 
runs as follows; we start from several "fairly large" graphs and construct all possible overlaps of their 
embeddings in the lattice in a recursive way. This limits the number of contributing graphs considerably, 
when compared to expansion (|J), but, except for regular lattices like the square or honeycomb lattice, 
still leaves the problem of determining g r and the contributions a r open. Indeed, for the square lattice 
where g r are rectangles, a r can be explicitly expressed via the ratio of the graph side lengths and the 
order to which the expansion is correct is in direct connection with the perimeter length of the largest 
graphs under consideration. For the quasiperiodic lattices which we wish to investigate the problem is 
not so simple. In what follows we focus on the PT Q and present the details of the expansion method 
for it in the next section. 



III. CALCULATION OF SERIES EXPANSION FOR THE PENROSE TILING 



The PT is an aperiodic tiling of a plane by two kinds of rhombi of unit length side with angles 2ir/5 
and 47r/5 respectively. A discussion of the methods of generation and geometrical properties of this tiling 
can be found in Jl3| , here we only mention a particularly useful feature, namely that embedding numbers 
of finite patches from this tiling can be calculated exactly and take the form n + mr where r = (\/2+ 1)/2 
is the golden number and n,m are rational numbers. The calculation of the series expansion consists 
therefore of the following steps: 
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1. Choose an initial set of "fairly large" graphs which are expected to be large enough that every 
connected subgraph of the underlying tiling with perimeter length not larger than a given threshold 
2L can be embedded in one of them. While on the square lattice this condition is satisfied by all 
possible rectangles with perimeter length 2L on the PT the things are worse due to the lack of 
periodicity of the tiling. Moreover, as opposed to the square lattice graphs in the PT can have 
different "boundary line fillings", i.e. there are different graphs having the same boundary line p3| . 
Knowing that the PT contains eight different vertex types, i.e. different site environments related to 
the nearest neighbours, we cut out appropriately large patches around each vertex type, obtaining 
eight patches, and found all possible "boundary line fillings" . There is still a lot of ambiguity in this 
procedure since a patch in not uniquely determined by the vertex type of its central site. It would 
be more correct to take all possible higher order vertex types p^ |, i.e. m-order vertex types related 
to neighbours located not further than m edges lengths from the site, but since their number grows 
quite rapidly with the order, the initial set of graphs would be too numerous and the generation of 
overlaps (see next item) too time consuming. 

2. Generate all possible overlaps of embeddings of the initial graphs in the tiling. 

It is difficult to estimate how the time of the generation depends on the number of initial graphs. 
Let us say a couple of words about this, however. We group graphs into generations so that the 
initial set of graphs constitutes the zeroth generation and the nth generation consists of overlaps 
of graphs from the (n-l)st and zeroth generations. Since the time for creating the nth generation 
depends on the product of numbers of graphs from generation zero, #g(0) and generation (n — 1), 
=tf=g(n — 1), starting from a too numerous zeroth generation should be avoided. The total number 
of overlaps grows rather slowly with #g(0) for large #<?(0) and most of the computing time will 
be devoted to checking and rejecting graphs which occurred before. On the other hand if we took 
too few initial patches, the covering of the lattice with them would be incomplete, there would be 
plenty of "holes" not covered by any of the patches, and thus the series expansion would be error 
laden. The rule of thumb is to take <?(0) not larger than twenty and choose the patches in such a 
way that their interiors differ as much as possible. 

Again, on the square lattice it's immediately clear that the overlaps are rectangles, because every 
rectangle can be constructed as an overlap of two other rectangles, whereas on the PT the shapes 
of graphs and their quantity depends on the initial set. To make the things worse we are not even 
sure that we obtain star graphs (page 1-16 in |l9|), i.e. graphs without articulation points, because 
the initial graphs are not necessarily convex. Connected graphs consisting of multiple components 
will cause some difficulties by the calculation of partition functions by the transfer matrix method 
(see following items). 

3. Calculate the contribution a r of graph g r in (0) in the following recursive way 



where the sum on the right-hand side runs over all graphs g v in which g r can be embedded. 
4. Calculate logarithms of partition functions log Z(g r ) by the transfer-matrix method. 



Here we have to distinguish two cases, namely the case when the graph has no articulation points (star 
graph) and the contrary (multicomponent graph). The latter is undoubtedly more complicated but 
fortunately it turns out that it takes place only in a minority of the graphs under consideration. Let us 
firstly discuss the case of a star graph. We can define a perimeter of the graph, i.e. a line consisting 
of edges each of which belongs only to one rhombi. The sum over spin configurations can be performed 
by moving a boundary line across the graph. At each stage the boundary line goes through a number, 
say k, sites. For the Q-state Potts model we have Q k different spin configurations on the boundary line. 
Now we define a Q k dimensional vector Z(cr) consisting of partition functions calculated for the patch 
composed of sites from the boundary line, with a given spin configuration er = {<J\, 02, . . . , 0jv} assigned 
to them, and sites already traversed by the boundary line. The initial values of Z(a) are given by: 




(4) 



The transfer-matrix method 



Z{a) 



i a (l-y) b where x = exp(-/3J), y = 1 - exp(-/3B/2) 



(5) 



and 



N-l N 

a=J2A(a p ,a p+1 ), 6 = 2^A(ct p ,0) 

p— 1 p— 1 



(6) 



Shifting the boundary line corresponds to generating a new vector Z'{a') of partition functions from 
the old vector Z . There is a lot of ambiguity in shifting the boundary line by a given number of tiles. 
In our case it amounts, however, to considering only three kinds of movements, by one tile, by two tiles 
and a shift between two given boundary line configurations, which we discuss in the following. Placing 
the initial boundary line on the perimeter of the graph and moving it at each stage by certain number 
of tiles, see FIG.0, we have performed the sum over all configurations after reaching the final position of 
the boundary line (also lying on the perimeter). 




FIG. 1. Shifting of a boundary line through a graph, from the initial (leftmost picture) to the final (rightmost 
picture), corresponding to calculating the partition function by transfer-matrix method. 

Now we discuss the details of updating the partition functions for the two kinds of boundary line 
movements, see FIG.||. 

1. One-tile movement 

For 1 < J < N we have 

Z'(a') = i a (l-y) b J2 (J) 

<TJ=l 

where a = A(a'j_ v a'j) + A(a'j,a' J+1 ) and b = fA(a'j, 0), / = 2. 

2. Two-tiles movement 

For l<L<P<Nwe have 

Z\a' p ) = i a (l-y) b ]T J2 ( 8 ) 

tTp + i — 1 CTp + 2— 1 

where a = A(a' L , <r' P+l ) + A(a' P+1 , cj' p+2 ) + A{a' P+2 , <j' L+l ) + A{a' Pl a' P+3 ) and b = fiA(a' P+1 ,0) + 
f2A(a' P+2 , 0) where fx = fa = 2 and the new spin configuration is permuted with respect to the old one 

a ' P = Wp^ a 'p^---^' PN ) and 




FIG. 2. Two kinds of movements of the boundary line, by one tile (left) and by two tiles (right). 
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3. Shifting the boundary line to the final position 



In most cases it is possible to displace the boundary line from the initial to the final position by a 
sequence of the movements defined above. Sometimes, however, we arrive in a dead end because none of 
the movements can be done, see FIG|| In this case we have to shift the line directly to its final position 
by summing over all the spins which have not been taken into account yet. The formal prescription for 
updating Z{a) in this case reads: 



Z'(a Fl , . . . ,cr Fg ) = 



E 



£ a (l-y) b Z(a Bl ,...,a B9 ) 



(10) 



where a = Yf j=1 A(a ejA ,a ej 2 ) and b = 2 A(cr Sj , 0) and B = {1,16,15,14,13,18,19,20,9}, F = 

{1, 16, 15, 14, 13, 12, 11, 10, 9}, e = {(9, 10), (10, 11), (11, 12), (12, 13), (12, 9), (9, 18)} and s = {12, 11, 10} 
denote the current and the final boundary lines, the edges and the sites which were not taken into account 
yet respectively. 




FIG. 3. A boundary line (middle picture) which cannot be pushed forward by performing one of the movements 
discussed above. The initial and the final line configuration are shown at the left and at the right respectively. 



Can the transfer-matrix formalism (tmf) also be applied to the case of a multicomponent graph? The 
answer is affirmative because every connected graph can be dissected into its star graph components for 
which the tmf is applicable. Since, however, star graph components share certain sites at their boundaries, 
which we call in the following isolated sites, we have to calculate a whole set of partition functions with 
given spin values at isolated sites and combine them to get the partition function of the whole graph. In 
the following we assume the simplest case namely, that every isolated site is shared by exactly two star 
components. This was indeed the case by our overlap graphs. Let us explain the procedure for the case 
of a graph depicted in FIG.|[ The partition function Z can be build up from partition functions Za{&i), 
Zb{?\i 02, 03), Zc[oz) and Z B (o2) corresponding to star components A, B, C, and D with isolated spins 
01, cr 2 and cr 3 . 

Z = ^2 Z A (cri)Z B (<J 1 ,a2.,a 3 )Zc(o- 3 )Z D (ij2) (11) 




FIG. 4. A multicomponent graph consisting of four star components (left) and the perimeter of a star graph 
with two isolated sites on the initial boundary line (consisting of sites from 1 to 7) and three isolated sites on the 
final boundary line (sites from 7 to 12) (right). The isolated sites are marked with circles. 

Now, the problem consists in calculating partition functions for a star graph with specified spins at 
isolated sites located at the boundary. Assume that we have p isolated sites jk, k = 1, ..,p located at the 



initial boundary line and q isolated sites Ik, k = 1, .., q at the final boundary line respectively, see FIG.|. 
The calculation of Z(sj 1 , . . . , Sj , s; 15 . . . , si ) amounts to repeating the tmf Q p+q times and modifying 
the initial and the final partition function set by setting certain entries to zero. We replace the initial 
partition function set by 



n ttokiSjk) 

Lfe=l 



Z(o\, . . . , cat) forgiven s^, . . . , Sj (12) 



and the final partition function set is multiplied by [Ofc=i ^"'in s i k )} again for a given spin configuration 
si t ,...,si . Another slight modification which is required consists in setting the factors and f% 
entering in the exponent b in equations ([^) according to whether the site is isolated (one) or not (two). 

IV. SERIES EXPANSION OF THE FREE ENERGY, MAGNETISATION AND FIELD 
DERIVATIVES OF THE MAGNETISATION 

We have performed calculations for a set of graphs constructed in the following way. We cut off seven 
fairly round shaped patches from the PT so that the central sites of the patches had different vertex 
types and their perimeter lengths were not larger than 30 edge lengths. Then we enlarged the set of 
patches by all possible "boundary line fillings" obtaining in effect twelve patches, see FIG.g. In the next 
step we constructed all possible graphs contributing to the expansion in the recursive way described in 



section [II. Their number turned out to be 1004. This part of computations was rather tedious, up to two 
weeks for the second set on a SunOS machine, because in generating graph overlaps many graphs turned 
up repeatedly and had to be rejected. In the next step we generated "decorations" of graphs, i.e. we 
determined vertex types of all sites of the graph including those on its boundary. Since the graphs could 
have several decorations the number of graphs we have to deal with increased to 5737. Now we were ready 
to compute the coefficients a r entering in (^) which appeared to be different from zero only for a small 
fraction of all graphs, namely for 154 graphs. This is not a surprising result since on the square lattice the 
vast majority of rectangles used in the expansion yields zero coefficients as well ^lj . Fortunately, most of 
the relevant graphs here were star graphs so we could easily compute the free energies log Z(g r ) entering 
in (||) in the way described in section III . There were however some awkward multicomponent graphs for 
which partition function computations were more tedious. The series expansion is shown beneath. 

After reordering the expansion (^), i.e. collecting together terms with the same power of y, the free 
energy F(x,y) takes the form. 

oo 

F(x,y) = logZ(0) = F (£)+F 1 (i)y + F 2 (i)y 2 + ... = ^ F n (x)y n (13) 

n=0 

Quantities like the spontaneous magnetisation M(x), susceptibility x(%) an d field derivatives of the 
susceptibility X {&) = d n x{x)/dy n can be expressed as linear combinations of the polynomials F n (x). 

M(x) = dF(Z,y)/dB\ B=Q =F 1 (x) 
X (x) = d 2 F(x,y)/dB 2 \ B=0 =2F 2 (x)~F 1 (5 : ) 
X (1) (£) = d 3 F(S: 1 y)/dB :i \ B=0 ^6F 3 (S:)~6F 2 (S : )+F 1 (S : ) 
X {2 \£) = d i F(5;,y)/dB 4 \ B=0 =2m(x)-36F 3 (i) + UF 2 (5 : )- F^x) 

(14) 
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FIG. 5. Patches from the Penrose lattice used as input for the finite lattice method calculations. The columns 
contain all possible "boundary line fillings" of seven patches the central sites of which correspond to seven different 
vertex types. 



53 54 55 56 57 5s 59 5io 5n 

M(x) -1.0832 -0.1803 -1.3736 -3.656 -10.0258 -12.8017 -11.2619 -16.135 -72.9404 

X(x) 2.1026 0.3607 4.4234 16.0886 63.0069 110.5216 137.4134 241.4194 990.5163 

X (1) (z) -4.0136 -0.7214 -15.448 -73.4949 -404.9629 -941.0482 -1484.2237 -3081.2191 -12969.5517 

X (2) (z) -33.5551 -5.861 -94.3 -405.998 -2039.6444 -4479.9119 -6807.027 -13838.423 -58117.8824 



512 513 514 515 516 517 518 519 520 

M(x) -231.8001 -493.5029 -841.3831 -1651.9956 -4125.5155 -9628.4929 -18432.4073 -33133.1842 -70657.8814 

X(x) 3397.510 8427.160 17816.643 41630.163 110568.346 277554.684 621098.121 1362853.02 3231583.007 

X W (x) -49115 -139934 -350428 -944924 -2738982 -7510844 -19047918 -47924567 -125178572 

X (2) {x) -217789 -612298 -1512009 -4036132 -11635964 -31747510 -79992423 -200008381 -520387269 

TABLE I. Expansion coefficients of the magnetisation M (x), susceptibility x( x ) an d its field derivatives x' 1 ^ {x), 
X^ 2 { x ) obtained from the finite lattice method. 



V. VERIFICATION OF CORRECTNESS OF THE COMPUTED EXPANSION 

There is a duality relation connecting the low temperature expansion of the Ising model on the lattice 
Q to the high-temperature expansion on the dual lattice T>, which takes the following form: 

Z g (x,y) = exp/3(MJ + NB)Zg(x,y) = 2 N (cosh (3 J) M (cosh (3B) N Z v (w, h) (15) 

where the low temperature variables are x = exp{— 2(3J}, y = exp{— flB} and the high-temperature 
ones are w = tanh{/3J}, h = tanh{/3£?}. In the field free case h = the high-temperature expansion 
of Z%>(w,0) can be expressed by the square root of the determinant of a 2M x 2A1 complex matrix 
|0,E3[ , which for periodic lattices amounts to calculating a finite-dimensional determinant the dimension 
of which is of the order of the size of the unit cell. Therefore the free energy expansion in variable x can 
be calculated by taking logarithms of equation (|15|). 

F = lim ^\ogZ g {x) = log2-flog(l-™ 2 )+logZ 23 H (16) 

where q = limjy >oo 2M/N is the mean coordination number. The expansion of the last term on the 

right-hand side 
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log Z v (w) = Y,9nw n (17) 

n=3 

is obtained from Kac-Ward determinants for large enough PAs of the Penrose lattice, see |u| for de- 
tailed explanation. In TABL E pi] we show the expansion coefficients g n for successive PAs together with 
the coefficients of Fq(x) (see |l3| ) obtained by the finite lattice method (FLM). The data for the high- 
est approximants are quite close to these for the FLM; the relative discrepancies for n = 3, ... ,20 are equal 
-4.9%,-0.2%,-0.4%,-5.1%,-4.7%,-2.8%,6.8%,14.5%,l.% -1.6% -6.4% -9.3% -11.1% -9.4% -8.2%, 
—8.4%,— 11.3%,— 13.1% and in most cases do not exceed ten percent. In addition both data sets depend 
on n in a similar way. Indeed, assuming known values for the critical point x c — 0.434269 and the critical 
exponent a = 2 we define the sequence r n in the following way: 

r n = 9n/g n -i ~ l/ar c (l - (a + l)/n) (18) 

This sequence approaches zero r n — ► for large n, see pages 187-199 in fl9|| . If we now compare the 
sequences from both the PA coefficients and the FLM coefficients we see that the relative discrepancies 
except for n=5, 10, 12 are all smaller than ten percent as well. 
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TABLE II. Expansion coefficients for free energy (|17|) for periodic approximants m — 3, . . . , 9 of the dual 
Penrose lattice with M edges in the unit cell. Underneath the solid line coefficients obtained by the finite lattice 
expansion. The numbers g n approach the expansion coefficients for the dual Penrose lattice when m — > oo. 



The lowest coefficients of our expansions can be also calculated exactly by counting graphs on the 
dual Penrose lattice. Here we confine ourselves to the free energy and the magnetisation expansions. We 
compute their first four nonzero coefficients and show that they are indeed close to those from tables 
TABLE | and TABLE [| We start from the low temperature expansion 

x M ' 2 y N Zg(x,y) = Zg{x,y) = 5>„,m^V m (19) 

with h n _ m counting graphs, in general multicomponent graphs, from dual lattice consisting of m sites and 
n bonds on the perimeter. It is readily seen from figures (FIG|| FIG.fi], FIG.^ and FIG.^j) that for the 
Penrose lattice the non zero coefficients take following values: 

fc M = (7 - At)N h ia = (-8 + 5t)N (20) 

h 5> i = (10 - 6t)N h 5>2 = (-16 + 10r)JV 

h 6)1 = (-21 + 13t)N h 6 , 3 = (-8 + 5t)N 

h 6 , 2 = (22 - 13t)N + (7 - At)N{{7 - At)N - 1) 
hi, i = (13 - 8t)N h 7 ,4 = (-8 + 5t)N 

h 7 ^ = (79 - A8r)Nh 7 a = (-63 + 39t)N + (7 - 4r)(-8 + 5t)N 2 



Let us notice that to coefficients ft.6,2 and h-j^ contribute also disjoint, two-component graphs thus the 
coefficients are second degree polynomials in N . 
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/ V 



/ ^ N S / 

/ \ X X / / 



r 



1 



(5 - 3x)N (2 - x)N 

FIG. 6. Graphs from Penrose lattice contributing to the coefficient /i3,i and their embedding numbers expressed 
through r = (\/2 — l)/2. The dual graphs are constructed by connecting midpoints of rhombi abutting at bonds 
terminated by filled circles. 
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FIG. 7. The same as above corresponding to coefficients /u,i(three on the left) and /is,i(last on the right). 
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FIG. 8. Graphs contributing to the coefficients he,3(le&) and /i6,i (right). 
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FIG. 9. Graphs contributing to the the coefficient /i6,2. 
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FIG. 10. Graphs contributing to the coefficient hi t i. 



i n 



* \ /V\ * 



2 x 2 x 1 x 1 x 

(26 - 1 6x)N (26 - 1 6x)N (26 - 1 6x)N ^ (26 - 1 6x)N 

' / 7 

\/*0, *\/*0, * 

' / / / > / / / > / / / > / / 

x // \ \ s // \ \ x // \ \ x y \ \ 
\ )2-x- * \ )2-x- * \ )2-x- * \ ) r - ^ 

N - '(-21 + 13x)N N - '(-21 + 13x)N ^ '(-21 + 13x)N N - '(13 - 8x)N 

' / X / / X 

\ \ / / x \/ / 

x/ ^ v s 

5 x 5 x 

(18/5- 11/5x)N (18/5 - 11/5x)N 



FIG. 11. Graphs contributing to the coefficient hi t z . 
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(13-8t)N (-8 + 5t)N 

FIG. 12. Graphs contributing to the coefficients 7i7,i(left) and /i7,4(right). 



If we now insert the coefficients from (21) into the definitions of the magnetisation M(x) and the 
susceptibility x( x ) 



M{x) 



1 d\og[Z g (x,y)] 



N dB 
1 dM{x) 



B=0 



1_ d\og[Z g {x,y)\ 
N V dy 



(21) 



N dB 



B=0 



1 d dlog[Zg(x,y)] 



dy 



we obtain following expansions: 

F(x) = 0.5279a; 3 + 0.0902a; 4 + 0.4721a; 5 + 0.8262a; 6 + 1.58359a; 7 - 
M(x) = 1 - 1.0557a; 3 - 0.1803a; 4 - 1.3049a; 5 - 3.4164a; 6 - 9.25233a; 7 - 
X(x) = 2.1115a; 3 + 0.3607a; 4 + 4.0526a; 5 + 14.6099a; 6 + 55.6843a; 7 - 
X (1) (x) = 4.2229a; 3 + 0.7214a; 4 + 13.8761a; 5 + 64.6563a; 6 + 341.449a; 7 - 
which conform quite well to the values from TABLE g. Indeed, the relative differ ences between both sets 
of coefficients do not exceed ten percent in any case. 



0(x s ) 
0(x 8 ) 
0{x s ) 
0(x 8 



VI. ASYMPTOTIC ANALYSIS OF THE SERIES EXPANSIONS 



Now the problem consists in extracting critical exponents from the obtained expansions. The simplest 
approach, the ratio method, in which one examines the asymptotically linear dependence of ratios g n /g n -i 
( |l8| ) on 1/n and obtains x c and a from linear regression, is inapplicable in this case because of the slow 
convergence of series. Indeed, the residues r n ( |l8|) are much larger than those for the square lattice and 
alternate in sign, see figures FIG.[ll| and FIG!^4|, what makes the asymptotic analysis difficult. This 
approach requires knowledge of x c which is known from other works Q , p4| only with a limited accuracy. 
Applying the Pade method gives much more satisfactory results. Assuming that our thermodynamic 
functions F(x) behave in the vicinity of the critical point x c like F[x) ~ (1 — x/x c )~ a A(x) it is readily 
seen that functions Go (a;) and G\{x) behave asymptotically as follows 

G Q {x) = ±. Lg^M)/-f (log^)) ~ £±1 + 0(x - x c ) (22) 
dx \ dx I dx a 



G\{x) =(x- x c )-^- QogF(x)) ~ a + 0(x - x c ) 
dx 



(23) 



Constructing Pade approximants to Go (a;) and G±(x) and evaluating them at x — x c = 0.434269 usually 
yields a reasonable estimation of q, even if x c is only known with a moderate accuracy. Results of this 



analysis, shown in tables TABLE III and TABLE [V are fairly close to the square lattice exponents. The 



sequence of Pade approximants is quite stable, except for few cases marked by asterix. Whenever the 
square lattice results converged Penrose data did as well. 

On the other hand we can compute biased estimates of x c , assuming known values of critical exponents. 
Indeed, the appropriate poles of Pade approximants to the function [^(a;)] 1 /" should give rapidly conver- 
gent sequence of estimates of x c . These sequences for the magnetisation expansion are shown in TABLE 
M. In most cases the data do not deviate more than one percent from the exact values x c = 0.434269 
(Penrose lattice |Q) and x c — y/(2) — 1 (square lattice). 

We can therefore claim that the data supports the claim that the quasiperiodic Ising model under 
consideration belongs to the square lattice (Onsager) universality class. 
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FIG. 13. Plots of residues r n (|l^) as a function of 1 /n for the free energy (left) and the magnetisation (right) on 
the Penrose and the square lattice respectively. In both cases we took the critical exponents a = 2 and /3 = 1/8. 
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Approximant 
















[n,n-l] 






[n 


,n] 






[n,n+l] 






Penrose 


Square 


Penrose 


Square 


Penrose 


Square 




Go (a:) 


Gx(x) 


Go(x) 


Gi(x) 


G (x) 


Gi(x) 


Go (a;) 


Gi(x) 


Go (re) 


Gx(x) 


Go(x) 


Gi(x) 


8 


0.159* 


0.128 


0.125 


0.125 


0.144 


0.126 


0.125 


0.125 


0.126 


0.126 


0.125 


0.125 


9 


0.242* 


0.126 


0.125 


0.125 


0.125 


0.126 


0.125 


0.125 


0.126 


0.127 


0.125 


0.125 


10 


0.121 


0.127 


0.125 


0.125 


0.123 


0.126 


0.125 


0.125 


0.126 


0.126 


0.125 


0.125 


11 


0.121 


0.126 


0.125 


0.125 


0.125 


0.126 


0.125 


0.125 


0.126 


0.124 


0.125 


0.125 


12 


0.140 


0.128 


0.125 


0.125 


0.135 


0.114 


0.125 


0.125 


0.113 


0.120 


0.125 


0.125 


13 


0.139 


0.121 


0.125 


0.125 


0.115 


0.118 


0.125 


0.125 


0.113 


0.122 


0.125 


0.125 


14 


0.067* 


0.128 


0.125 


0.125 


0.042* 


0.097* 


0.125 


0.125 


0.131 


0.137 


0.125 


0.125 


15 


0.069* 


0.154* 


0.125 


0.125 


0.128 


0.146 


0.125 


0.125 


0.130 


0.100 


0.125 


0.125 



TABLE III. Estimates of the magnetisation critical exponent (3 by means of Pade approximants [n, m] to 
functions Go(x) and Gi(x) (|2^) constructed from the expansion to order 20 for the square- and Penrose lattice 
respectively. Entries marked by asterix differ strongly from square lattice exponents. 
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Approximant 
















[n,n 


-1] 






[n,n] 






[n,n+l] 






Penrose 


Square 


Penrose 


Sc[UcL 


re 


Penrose 


Square 




G (x) 


G 1 {x) 


G (x) 


Gi(x) 


Gq(x) 


G x {x) 


G (x) 


(j\\X) 


Go («) 


Gi(a;) 


Go(x) 


Gi(a;) 


8 


-1.287* 


-1.855 


-1.383* 


-1.743 


-1.469 


-1.778 


-1.741 


-1.740 


-1.095* 


-2.230 


-1.741 


-1.741 


9 


-1.313* 


-1.832 


-1.741 


-1.741 


-1.397* 


-1.946 


-1.741 


0.000* 


-0.808* 


-2.244 


-1.699 


-1.585 


10 


-1.299* 


-1.834 


-1.702 


-1.585 


-1.723 


-1.790 


-1.716 


-1.454 


-1.259 


-1.559 


-1.716 


-1.587 


11 


-1.431 


-1.851 


-1.716 


-1.587 


-1.133* 


-1.249* 


-1.716 


0.000* 


-1.088* 


-1.758 


-1.732 


-1.747 


12 


-1.095* 


-3.743* 


-1.738 


-1.747 


-1.131* 


-0.244* 


-1.736 


-1.747 


-1.259* 


-1.129* 


-1.736 


-1.747 


13 


-1.472 


-1.569 


-1.736 


-1.747 


-2.005 


-0.956* 


-1.736 


0.000* 


13.636* 


0.316* 


-1.757 


-1.748 


14 


-0.336* 


-0.752* 


-1.737 


-1.747 


-4.638* 


-2.848* 


-1.742 


-1.748 


-4.610* 


-2.432* 


-1.742 


-1.748 


15 


-4.610* 


-2.458* 


-1.742 


-1.748 


-4.637* 


-1.924 


-1.742 - 


-0.001* 


-7.190* 


-2.411* 


-1.743 


-1.749 



TABLE IV. As before for the susceptibility critical exponent 7. 
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[n,n-l] 






Approximant 
[n,n] 


[n,n- 


HI] 






Penrose 




Square 


Penrose 


Square 


Penrose 




Square 


8 


0.434 




0.414 


0.462 


0.414 


0.437 




0.414 


9 


0.433 




0.414 


0.434 


0.414 


0.434 




0.414 


10 


0.434 




0.414 


0.434 


0.414 


0.434 




0.414 


11 


0.434 




0.414 


0.434 


0.414 


0.435 




0.414 


12 


0.433 




0.414 


0.446 


0.414 


0.437 




0.414 


13 


0.434 




0.414 


0.435 


0.414 


0.436 




0.414 


14 


0.431 




0.414 


0.438 


0.414 


0.457 + 0.018 *7 




0.414 



TABLE V. Biased estimates of the critical point x c — exp(— 2/3 c ) obtained from the magnetisation expansion 
to order 20 for the square- and Penrose lattice respectively. 



VII. CONCLUDING REMARKS AND OUTLOOK 



The aim of this work was to analyse quasiperiodic Ising models by means of graphical series expansions. 
We calculated low temperature expansions of the free energy, magnetisation and the susceptibility to order 
20 and extracted the respective critical exponents. We note that we did not obtain exact values for the 
exponents. They are marked by errors, which however do not exceed 10% in most cases. This feature is 
rather unusual for series expansions on regular lattices where the coefficients are exact to a given order 
which is determined by the size of patches used for calculations. This deserves further comment. There 
are the following two sources of errors in our method: 

1. There are always graphs which cannot be embedded in any of the patches used in the calculations. 

2. The FLM consists in covering and probing the lattice with a finite set of patches. Since the lattice 
is irregular, in particular, not periodic there can always be "holes" , i.e. groups of sites, associated 
with graphs on the dual lattice, which are not covered by any of the patches. Such incomplete 
covering with only one single patch is shown in figure [l5| Consequently, even if the patches are 
large there can always be few small graphs which cannot be embedded into them and which produce 
an error in even the lowest order coefficients. 



1 1 




The next step in our research is to analyse the quasiperiodic Q-state Potts models, especially for Q = 
3, 4 because the Harris-Luck criterion implies that quasiperiodic order should be strong enough to alter 
the critical behaviour in these cases. The main problem here consists in calculating partition functions 
for finite patches. This can be done, for instance, by improving the FLM |^5| where the expansion for a 
particular Q is obtained from partition functions with smaller Q values. Another possibility is to use the 
Fortuin-Kasteleyn representation of the Potts model |17|, expressing partition functions as a function of 
Q. Work in this direction is in progress. 
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